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We critically review the fast algorithms for the numerical study of two-dimensional Josephson 
' junction arrays and develop the analogy of such systems with electrostatics. We extend these 

, procedures to arrays with bus-bars and defects in the form of missing bonds. The role of boundaries 

and of the guage choice in determing the Green's function of the system is clarified. The extension 
of the Green's function approach to other situations is also discussed. 
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I. INTRODUCTION 



The dynamical properties of Josephson junction arrays (JJAs) are currently the foci of several experimental and 
theoretical investigations These arrays can now be routinely fabricated in several sizes and geometries, and the 
characteristics of their junctions can be varied at will over a wide range of values 0]. A large body of high-precision 
experimental data has consequently become available for JJAs in the presence of external magnetic fields |3|-|5| ■ On 
the theoretical front, several insights into the behaviour of JJAs have come from numerical studies of the underlying 
equations of motion as given by the resistively- and capacitively-shunted junction (RCSJ) model using input current 
drives and defects, both controlled Q and random [Q. With the size of experimental arrays increasing continuously, 
and with the number of interesting effects best seen only in large arrays going up in equal measure, it has become 
imperative to find ever more efficient algorithms for implementing the corresponding simulations, inclusive of all the 
experimental conditions. An example of the latter for current-driven arrays is the presence of bus-bars, through which 
. the external current can be conveniently injected or withdrawn. 

To understand the problem which these algorithms must address, we recall that in the RCSJ model, the total current, 
iij, (inclusive of external drives where applicable) flowing through the junction between sites i and j, is viewed as 
O ■ consisting of three 'channels' in parallel: superconductive, resistive ( or ohmic) and capacitive. The currents in each 
^ ', of these channels can be expressed in terms of the phase difference, 9ij = &i — 9j, across the junction. This leads to 
? the following equation for the evolution of the latter in time: 

Here C, R and I c are the shunt-capacitance, shunt-resistance and critical current of the junction respectively. We 
note that Eq. (|l|) holds under assumptions of zero temperature, zero magnetic field and infinite perpendicular magnetic 
penetration depth. The last of these allows us to neglect self-induced magnetic fields. 

Using total current conservation (TCC) at each site ||, we can write 

5>^ + ^ +sin%= I- Vi (2) 
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where If xt = if xt /I c is the normalised current being fed to or extracted from the array site i (see Fig. ), j3 c = 
2eI c R 2 C/h is the McCumber-Stewart parameter and r = t(2eI c R)/h is the time measured in units of the characteristic 
period lo~ x = Ti/(2eI c R). The summation on j, over the nearest neighbours of i, can alternatively be expressed in 
terms of a multiplication by the discrete laplacian, Gq 1 . Eq.Q then assumes the matrix form 



(Gq 1 ) . h = ^[ir - V kj + sin%] = -di (3) 



J {ij) 



If we set 9i — Xi and 9i — Ui, the complete set of dynamical equations reads (Go = — ^(Mi [?/]) an d = Z/i- 

For the overdamped case, corresponding to f3 c = 0, the relevant equations are 

E( G o = (4) 

j 

where — di is now given by If xt — sinfljj but Gq 1 is, of course, the same as in Eq.(H). 

It follows that for an N x x N y = N array each integration time step of Eq.(|]) has a complexity 0(N 2 ). This is 
because at every upgradation of the TV state variables, 9, in the under-, and 9, in the over-damped case, the constant 
N x N matrix Go has to be multiplied by the divergence vector [d] . It was first noticed by Eikmans and Himbergen 
H that the form of the Gq~ 1 is such that this multiplication can actually be carried out in 0(N In N) steps. The 
procedure was subsequently improved upon by Dominguez et al. who combined the fast-fourier transformation 
used by Eikmans et al. with the method of cyclic reduction to achieve a roughly 30% increase in speed. 

It is noteworthy that these algorithms are applicable even in the presence of an external magnetic field. Indeed, the 
application of a field, Bqz, perpendicular to the array transforms the phases 0y into the gauge-invariant combinations 
9ij + 2tt Aij where Aij — 1 /4>o y~]i A - dl, A is the vector potential and 4>o = h/(2e) is the elementary flux quantum 
threading a plaquette. This transformation clearly affects only the divergence term, and leaves untouched the matrix 
Gq" 1 , on whose form the algorithms are based. Similarly, white noise, which is taken into account by introducing 
a noise current into Eq.([|), also modifies only the divergence and hence does not affect the applicability of these 
algorithms. The statement continues to hold even if we grade each bond, i.e. make the R and I c junction-dependent 

In this paper we extend these fast algorithms firstly to arrays with busbars and secondly to those with defects in the 
form of missing bonds. As has already been mentioned, busbars often form the current injection and/or extraction 
edges of experimental arrays. In experiments involving vortices, which are repelled by busbars, these have been 
used to produce collimated vortex-streets (TTJ . The dynamics of vortices have also been investigated with one edge 
shorted by a single busbar fl2|| . Arrays with defects have likewise arisen in a number of contexts. The breakdown of 
superconductive flow in current-driven arrays with linear defects, for example, has been investigated in some detail. 
The exploration of the multiple- vortex sector, which arises in this study, requires running on large arrays and is all 
but impossible without the algorithm we discuss below p3| . Defects can also be used to provide a collimated beam 
of vortices . 

The paper is organised as follows. In section H, we present a simplified derivation of the Eikmans-Himbergen 
algorithm and interpret all the key equations in the familiar language of electrodynamics. Apart from being more 
transparent, our derivation clarifies the role of boundaries and the connection between the lattice and continum 



descriptions of the systems being studied. These insights are used in section III to generalize this algorithm to arrays 



with busbars. The case of single busbar forces us to resort to the technique of cyclic reduction, which we consequently 



discuss at this point. In section IV, we extend these algorithms to JJAs with defects, created by eliminating or adding 
bonds. Section ^ contains a summary and discussion of our results. 



II. THE EIKMANS-HIMBERGEN ALGORITHM 
A. Preliminaries 



All efficient algorithms for the numerical integration of Eq. Q make essential use of the properties of the discrete 
laplacian, Gq" 1 . We thus begin by listing the more important properties of this matrix. From its definition (as 
given by Eqs. <^ and ^), it follows that Gq 1 specifies the connectivity between different sites. More precisely, 
(Gq 1 ) i . = (Gq 1 ) j { = —I for i 7^ j (site i connected to site j) and otherwise. The diagonal element (Gq = total 

number of points to which the site i is connected, and hence Tr(G ~ 1 ) is twice the number of bonds in the network. 
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Clearly, Gq 1 is real and symmetric. As a result, its eigenvalues are real and its eigenbasis complete in the N- 
dimensional space of states. Furthermore, det(Gp ) = 0. This can be deduced as follows. We first note that 
since Eq.(j^) involves only the phase differences (and their time derivatives) it is invariant with respect to the global 
transformation 9\ — > 9\ + cx(t) Vi. Substituting this into Eq.(|J), we have 

E ( G o%& + aw) = -* = E «■•..'>./', 

This implies that c\{t) yV (Gq~ ) „ = Vi. I.e., if we fix i and sum over all the columns j, we get zero. (The same 

is true, of course, if we fix j and sum over i since the matrix is symmetric). Hence det(Go~ 1 ) = and only (N — 1) of 
the N equations in Eq.(|J) are independent. (We could, alternatively, note that ('l'o)i = lVz is an eigenvector of Gq 1 
with eigenvalue 0). 

To explicitly evaluate 9i from 9i we have to eliminate the extra degree of freedom. For what immediately follows, 
the most convenient choice is Yli @i — 0- Eq.(Q) is then replaced by 

E(^o~%^ = -A (5) 



j 



where A = i = 1 • • • (N - 1) and D N = while (Gq 1 ),^ = (G^ 1 )^ i = l---(iV-l), j = 1---JV and 

(^o" 1 )^ = IVj- Moreover, the vector [d] must satisfy di = as can be seen by performing an additional sum over 

i in Eq.(^) and using the fact that J2i {@o ~ 0- For JJAs this condition is automatically satisfied since 9%j and 
sin 6ij are odd functions of their arguments and the net external current fed to the array is zero. 

One can, in principle, invert Qq 1 and determine 9 from [D\. This is the usual 0(N 2 ) process. We shall refrain 

from adopting this procedure and return, instead, to working with Eq.(^|). 

In doing so, we shall find it useful to keep in mind some electrostatic analogs of the equations we happen to be 
dealing with. These emerge clearly if we write the current conservation equation at the site k as 

I^=Il,+Il,+lC (6) 

Here fi=l,2 represents the x and y directions respectively. Furthermore, I^ff" 1 , ^k,a = s ^ n ^^k, I^u = and 
Iji? 1 * are the total, superconducting, normal and external currents, respectively, at node k in the p th direction. Using 
the TCC condition, V • I\ otal = 0, we have 



(Tl + T? 1 ) (7) 



On comparing Eq.(^) with the Laplace Equation V 2 (j> — —ptotai/^a and recalling that ptotai — Pfree + Pbound, we 
arrive at the following correspondences 

D(r) = ir l (8) 

-P(r) ee ff (9) 

<f>(r) = 9\ (10) 

Thus the condition imposed above on the 9\ is nothing other than a choice of reference potential or gauge. We, 
furthermore, note that the JJA thought of as a dielectric medium is highly non-linear. Indeed, Pk,n = sin(— J Ek tfl dt), 

where Ek 4l — A^9k is the p th component of the electric field at site k. This is to be contrasted with P(r) cx E(r) for 
a linear dielectric. 

B. The O(NlnN) procedure 
Going back to Eq. (|3j) , Gq 1 can alternatively be inverted by defining the Green's function [see Jackson] 

G(r,r')= E (11) 
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where Ag are the eigenvalues and ^(r) the orthonormalised eigenvectors of G , i. e. 

<% 1 n(?)=w?) (i 2 ) 

and r are points on the lattice. More explicitly, the site coordinates (xt, yi) = fi referred to an origin located at the 
lower left hand corner of the lattice, are related to the site indices i used earlier by i = Xi + yiN x + 1. We note that 
to get a well-defined function, we must explicitly remove from the sum in Eq. ([□]), the zero mode, which necessarily 
exists, since Gq 1 is singular. 

The removal of this mode, however, creates a slight problem. Indeed, if we premultiply by G^ 1 and use Eq.(|l2|) we 
get 



Y,G,\r,r")G{r"^ = £ ^A^r^r*') (13) 
f" k^O k 



■\/ L X Ly l ' J yj L X Ly 



which in discretised form is 



Now if we multiply Eq.(Q) from the left by G we get d\ — 'Y i L ' which reads just right, viz. 

6i = G l3 d J (15) 

provided the Q\ satisfy 9i = 0. In other words, with an appropriate gauge choice, G does allow us to invert Gq 1 . 
Indeed, the requirement that it does so fixes the gauge uniquely. It is to be noted that Eq.([l5|) involves [d] and not 
[D]. 

The function, G(r, r '), can be easily evaluated for a periodic lattice. The form of the operator G l (f, f '), in this 
case, is given by 

Go 1 ^, r') = A5 PjPI - 8?±£ xt f< - 6?± Syt ?< (16) 
It is easily checked that the normalised eigenvectors are of the form ex P lk ' r and correspond to the eigenvalues 

Ajj = 4 — 2 cos k x — 2 cos k y (17) 

We note that Xj: =0 = 0. 

Due to the rectangular periodicity of the lattice the eigenvectors must satisfy ^dr) — ^t(f+ LiBi) i = x,y. To 
fulfil this condition, the wavevectors must in turn be quantised as ki = (2mr)/Li i = x,y and n = 0, 1, • • • (L;_i). 
The Green's function of interest is consequently 

G(r,r') = — !— V " expifc-(r-f') (18) 

K ' L X L V ^ 4 — 2 cos k x — 2 cos k v v K ' y ' 

We can now estimate the number of steps required to evaluate Ylf G(f,f')d(f'). The sum over f is a fourier 
transform and can be carried out in NlnN steps to produce d(k). Then come the N multiplications X^d(k) followed 
by the sum over k, which has the form of an inverse fourier transform, and requires an additional NlnN steps, making 
a grand total of N(2lnN + 1) multiplications. 

Quite generally then, 6(f) can be evaluated from Eq.(|l5"|) by (i) creating an L x x L y matrix Go(k) = Nz/A? where N^ 
are the normalisation factors due to ^e(r), (ii) evaluating d(k) = J2r' d(f')^ , this being the forward transform 
(W), (iii) computing d'(k) — Go(k)d(k), and (iv) taking the 2-D inverse transform W to get 9(f) — $ (k)^> Af) . For 
certain types of transforms and certain values of L x and L y , W and W can each be performed by matrix decomposition 
or row-column techniques with a complexity 0(N In N). We now turn to a discussion of this procedure in the context 
of finite arrays. 
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C. The Imposition of Boundaries 



For a finite array, the gener al form of Gq 1 (Eq.©) is modified since now fewer than four bonds meet at all sites 
on the boundary (see Fig.(l)). For non-corner points on the left edge of the array for example, the operator has the 
form 3<5f?f — ^f+g^y — &?±e y ,r'- However, if \& satisfies f (xo,j/) — ^(a^o — 1, i/) = Vy along this edge, one can 
continue using the periodic form of G^ 1 because then (iSpf — S?±g x ,F> — Sf±g fi)^{xo,y) automatically reduces to 



conditions on if? 



5p±e t fti)^(xo,y). We thus see that the finiteness of the array imposes the following boundary 



*(xo,y)-*(xo-l,y) = Vy 

y(x,y )-y(x,y o -l)=0 Vx 

®(x L ,y) - ^(x L + l,y) = Vy where x L = x Q + L x - 1 

*(x, - *(x, y L + 1) = Vr where y L = yo + £ B - 1 



(19a) 
(19b) 
(19c) 
(19d) 



where (xo,yo) and (xl, Vl ) ar e the diagonally opposite corners of the rectangular lattice. It is amusing to note that 
the left-hand side of Eqs.(19a- 19d) are the discrete derivatives of ^f?(x,y) along the edges x — x n ,y — y n ,x — xl and 
y = y L , respectively. Thus the finite case poses for us a discretised Neumann boundary value problem. 



: pot 

Furthermore, we note that Ag (see Eq.(|17j)) has the symmetries of the square, i.e. A 
Thus any linear combination of the corresponding eigenvectors, 

^(x, y) = a\ expi(k x x + k y y) + a,2 expi(—k x x + k y y) 
+03 exp i(k x x — k y y) + 04 exp i{— k x x — k y y) 

is an eigenvector of Gq 1 with the same eigenvalue Ag. To find the specific linear combination, which satisfies the 
boundary conditions given above, we need merely impose each of Eqs.(19a-19d) in turn. Eq.( |l9a| ) can be satisfied by 
choosing xq = 1/2, a\ = 02, 03 = 04 and the resulting wavefunction is 



A 



A_fc_._t 



(20) 



if?(x,y) = 2a\ cos k x x exp(ik y y) + 2a,3 cos k x x exp(— ik y y) 



(21) 



Subjecting Eq.(|2l|) next to the constraint Eq.(19t), we get 

*f?(x, y) = ia% cos k x x cos k y y (22) 

where yo is now required to be 1/2 as well. The fact that xq = yo = 1/2 means that the origin of coordinates is chosen 
on the dual lattice and both coo rdina tes of eve ry la ttice point are half integers. 

The boundary conditions Eq.(19c) and Eq.( 19d ) at the right and upper edges of the array can be satisfied by 
imposing a quantisation condition on k. Indeed, using Eq.(19c) we have 



3ai cosfc y ysin(fc 2; i 2 ;) sin( — ) = Vy 



(23) 



whereby k x = n x n/ L x , n x = 0, 1, • • • , L x — 1. Similarly using Eq.( 19d ) k y = n y Tt/L y , 
Thus finally, the orthonormalised wavefunctions we are looking for are 



0,1, 



L x L y 



1 - 2^*.° 



1 



h y ,o cos(k x x) cos(k y y) k ^ 



The resulting Green's function to be used in Eq.(ll) is thus (^] 



G(r,r') 



E 



(1 



J k x ,0 



L x L y rr* 1 4 — 2 cos k x — 2 cos k y 



cos(k x x) cos(k y y) cos(k x x') cos(k y y') 



(24) 



(25) 



k^O 



The action of this Green's function on a uniform current drive, I ext (x,y) = I ext (S x 
shown to be jl5| 

G(x, y, x', y')I{x', y') = J(^y - x) 

y' 



x a - 8 x ,x Lx ) can be rigorously 

(26) 
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III. INCLUSION OF BUS-BARS 



Using the techniques of the previous section we generalise the procedure to arrays with either a single or a double 
busbar (henceforth referred to as SBB and DBB respectively). 

In the SBB case, the busbar is normally placed along the current-extraction-edge (see Fig. (2)). The corresponding 
9 and 9 are zero for all time i.e. they do not evolve. The drive edge is however open and can be connected to any 
kind of current profile. Such an array, with a linear profile, has been used previously to study vortex dynamics [ fL2"| . 
On the other hand, an DBB array has shorts along a pair of parallel edges (say those at x = xo and x = xl + 1 as in 
Fig(3)). In this case, no current can be injected in the x- direction but the y— direction is available to an arbitrary 
current drive. Such an array with electrically connected busbars could be used to simulate e.g. the "channel" in the 
experiments conducted by Van-der-Zant et al. [0. 

We note that by setting the 9 along the busbar to zero and measuring all other 9 s with respect to these preassigned 
variables, we are unambigiously fixing the gauge. Furthermore, since the wavefunction in a rectangular system is 
separable, i.e. ^(x,y) = 'I/x (x)\E f 2 (y) , any change in the boundary conditions along the x— direction, has no impact 
on ^(y), which therefore continues to have the form derived in Setion II, viz. 




*2(y) = \ ( 1 - ^kyfij cosk y y 

The altered b.c.s however enter crucially into the determination of &i(x) and hence in to that of ty(x,y). 

A. Arrays with a Single Bus-bar 

For the SBB case, we have the following b.c.s 

*i(xo) - *i(z - 1) = Vy (27a) 



^i(xl + 1) = Vy where xl — xq + L x — 1 



(27b) 



The second of th ese results from the fact that the wavefunction must vanish at the shorted edge as d iscus sed above. 
Condition Eq.(27a) fixes the form of ^i(x) to be ~ cosk x x where a; is a half-odd integer, while Eq.(27b) quantises 
k x - k x — ((2n x + / (2L X + 1) n x = 0, 1, • • • , L x — 1. The normalisation constant is fixed by the requirement that 
A 2 Ep=i^s 2 k x (p-^)^l. Thus 



*i(a:) 



- I cos k x a 



(28) 



and hence the Green's function G for the case of a single busbar is 



G(f, f" 



(1 - X,o) 



L y {L x + |) 4 — 2 008^ — 2cosfc y 



cos(fc.j;a;) cos(k y y) cos(k x x') cos(k y y') 



(29) 



This G coincides with the Go for this case because the zero-mode is absent. 

The action of G on a uniform input drive I ext (x, y) = I ext 5 x Xo can, as in the case of Eq. (|26|) , be exactly shown to 
be Hi 



G{x, y, x', y')I(x', y') = I(L X -x+\\ 



(30) 



Although a ID Fast Cosine transform can be used along the y— axis, this cannot be done for the ^-direction, since 
no fast Cosine transform of the form 



L-l 



C{m) = C( n ) cos 



n=0 



(2n + l)(2m+ 1)tt 
(4L + 2) 



G 



is known to exist |16|. If we resort to usual matrix multiplication along e x , the algorithm becomes 0(2NN X + 
2N In N y + N) in complexity and is useful only if N y >> N x . 

A much more efficient approach is to combine the FCT along e y with cylic reduction along the x-direction. This 
can be carried out as follows [0. We begin by writing Eq.(^) as 

@x — l.y ~t~ 9 xy tT y ' y -+- Qx + l.y — d x , y (31) 

where the operator T v > y is defined to be T y i tV — 8 y >, y ~i + S y >, y +i — 45 y > , y Taking the Cosine transform along the y- 
direction i. e. setting 6 x _ y = ^2 k Q x .k y cos k y y, we find that 

Sx,y'T y/yV = ^ 2 ( cos k v - 2 )®x,k y cos k y y (32) 

ky 

As a consequence the y- cosine transform of 9 X iV satisfies 

9X-I,ky + Hky)0x,ky + 0X + I,ky = ~ ^X.Jfcj, (33) 



where we have used X(k y ) = 2(cosfcj, — 2). Writing Eq.(|33J) thrice with x set equal to x — 1, x and x + 1 respectively, 
multiplying the second of these equations by — X(k y ) and adding all the three of them together, we get 

0x-2,k y + X w (k y )9 X:ky + 6 x+2 ,ky = -d x -i t k y + X(k y )d Xtky - d x+1<ky = -d^ ky (34) 

where X^(k y ) = 2 — (X(k y )) 2 . We note that the x values occurring in Eq.(|4|) increase in steps of 2 as opposed to I 
in Eq.(33). The two equations are, however, of the same form. We can thus use the reduction procedure repeatedly 
to get equations with x-values increasing in steps of 4, 8, • • • until after m stages (2 m = L x ), we end up with a single 
equation for the central line of variables: 

X^ l) (ky)9 Xa+L:i;/ 2.ky = di"+ii/2,fc y - 0X O ,ky - dx L +l,ky ^ 

The X^(k y ) and d { ^ k occurring in these equation are defined by the recursive relations: 

A (p) (M = 2-(A^(fc y )) 2 (36) 

41 = " A(P_1) 3 + d^l Kky (37) 

Thus, if XQi k y and XL+ \^ ky are known, then 6 Xo -\-L x i2,k can be immediately deduced. The knowledge of Q Xa +L m ii,ky 
leads to that of 6 Xo+Lx /^ ky and 8 XQ+SLx /^ ky which in turn determine Xo+qLx / Stky = 1,3,5,7) etc. As for the 
starting values, Q XL +\ lky — V k yi since this is the y-cosine transform of x , y for x = xl + lj i-e- taken along 
the busbar. Xo , ky must however be determined by explicit multiplication. Fortunately, since we know the relevant 
Green's function explicitly, we can carry this out in 0{N) steps, by taking the cosine transform of Eq.(|l5|): 



y x ,k v 



-TT , 2 — — cos k x xn cos k x x'd x i k (38) 



L y (L x + i) ^ *—f 4 — 2 cos k x — 2 cos k y 



(p) 

The number of steps required to evaluate all the higher-order divergences a£ k , p = 1, • • • m— 1 used by the method 



can be shown to be N y {(N x /2) — hi2 N x } 15 . With these divergences in hand, the cyclic determination of 9 x ,k y 
has a complexity O(N). Adding to all this, the number of multiplications involved in carrying out the forward and 
backward FCTs in the y-direction, we conclude that the entire procedure can be carried out in 0(2NlnN y + 2N + 
N/2 - N y ln 2 N x ) steps. 
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B. The Case of a Double Bus-bar 



Turning now to the DBB case, we note that the wave- function, must now satisfy the following b.c.s 

*i(a;o)=0 Vzj (39a) 



¥i(x£ + l) = Vy 



(39b) 



Condition Eq.(39a) fixes the form of ^i(x) to be ~ sink x x where x is an integer while condition Eq. ( p9b| ) restricts 
k x to the values (n x Tr)/(L x + 1), n x = 1, • • ■ L x . The normalisation constant, A, is fixed by A 2 JZp=i sin 2 k x p = 1. 
The resulting orthonormalised wavefunction is then 



and hence G for an array with a double bus-bar is 



L x + 1 



G{r,r') 



L y (L x + 1) 4 — 2 cos k x — 2 cos k y 



sin fc^x 



sin^x) cos(kyy) siri(k x x') cos(k y y') 



(40) 



(41) 



In this case, the transformations along both directions are amenable to fast, i.e. O(NlnN), procedures provided 
L x + 1 and L y are both integral powers of 2. 



IV. DEFECTS 

We move on next onto arrays containing defects in the form of broken bonds. As mentioned in Sec. |], arrays of 
this sort have recently attracted the attention of several groups. Xia and Leath showed that in contrast to networks 
of resistors p8| , where breakdown emanates from the most critical defect, J J As with missing bonds do not become 
resistive, the moment any one junction turns critical. Cohn et al. p9]| showed that defects exert an additional pinning 
force on a vortex placed in the system. Datta et al. [ p"3| demonstrated the existence of multiple vortex sectors in the 
steady state configurations of the system. They did so by running on large arrays using a modified version of the 
above algorithm which we now describe in detail. 

The evolution equation 9 for an array containing a single missing bond between the sites k = (x , yo) and (k + 1) = 
(xq + 1, yo) can be written in matrix notation ( see Eq.(g)) as 

{G^ + h)[6] = [d\ (42) 

where Gq 1 is the discrete Laplacian with free, periodic, SBB or DBB boundary conditions and h is given by 

hu.k = fofc+i,fc+i = 1 (43) 
hk,k+i = hk+i,k = — 1 

hij = i, j =/= k, k + 1 

Multiplying Eq.(^2|) with G from the left we get 

(I + A)[0] = G[d] = £ (44) 

where A = Gh and I is the identity matrix. The process of determining G[d] — £, as already been outlined and can 
be executed in N\nN steps. One then needs to evaluate (1 + A)~ £ efficiently. To this end, we note that A has the 
form 

Aij = Xi(6j,k - Sj,k+i) (45) 
where the numbers xi are got by explicitly by multiplying Go and h. From Eq.(^5|) it follows that 

A q = A(Ae) q ^ (46) 
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where Ae = Xk — Xk+i- In fact, Ae is the only non-zero eigen value of the matrix A and for the case discussed in 
Section |ll C| it is given by 



Ae = -4- V f lSk f f kv ;y (cos(fc x (x + 1)) - ^(fc^o)) 2 cos 2 (fc,y ) (47) 

Lt rLtqi 4 — Z COS Kx — ^ COS rZu 

where hi = (riiir)/ Li, m = 0. . .Li — 1 and i — x,y. From Eq.([l6|) it follows that 

(/ + AY 1 = I- A + A 2 -A 3 -\ = I ^- (48) 

v y 1 + Ae v ; 

Since the constant matrix R = —A/(l + Ae) contains only two non-zero columns, which are moreover negatives of 
each other, [9] can be determined from £ in just N multiplicative steps. 

The extension of the above procedure to the case of a linear defect Q consisting of n broken bonds is straight 
forward as each additional defect introduces into h a 2 x 2 diagonal block of the form given in Eq.(|4^). In this case, 
the series (Eq.(^8|)) cannot be summed analytically due to the presence of cross terms. An analytic summation is, 
however, unnecessary. It is sufficient to note that 

A i:j = x\{5 ]M - 5 jM+ x) + xj(S jjk2 - 8 jM+1 ) + ■■■ xV-(Sj, kn - S^ kn+ i) (49) 

for bonds missing between ki and fc^+i (i = 1, 2 • • - n), and that consequently A q has the same form as above. The 
numbers x^'s occurring in A q are of course g-dependent. It follows from this that the series can once again be written 
in the form 

(I + A)~ 1 =L + R (50) 

where R is a constant matrix of the form given in Eq. ([49]) . Clearly the action of R on a vector of length N can be 
determined in nN multiplicative steps. The overall complexity of the algorithm in this case is thus 0(N In N + nN). 
Since the number of broken bonds required to observe various breakdown phenomena |6| is much smaller than N, this 
results in a very large saving in the time required to perform the computations. 

It is interesting to note that the Green's function for the broken bond system G bb satisfies the Dyson's equation 
@ 

G bb = G + G(-h)G bb (51) 

The matrix (— h) is thus the analog of the potential due to an impurity in an otherwise perfect lattice and the case 
of the single defect is the analog of the S function impurity at a given site in the system. 

The procedure developed above is for the situation where no pair of broken bonds has any site in common. This 
ensures that the block matrices introduced in h are always 2x2 rather than m x m, m > 2 The general case can be 
dealt with in a similar manner albeit with an increase in complexity. Lastly, bonds can be added in the interior of 
the network rather than being eliminated since this does not violate the constraint dt = 0. 



V. SUMMARY AND DISCUSSION 



To summarize, we have extended the Eikman-Himbergen algorithm to a number of practically occurring situations. 
In particular, we have extended it to the case of a bus-bar placed along one edge of a rectangular array as also to 
that of electrically connected bus-bars placed along two parallel edges. Finally we have shown that our simulation 
can be carried out, with only a marginal increase in complexity, for bonds eliminated from or added to the array. In 
all cases, the open edges can be connected to current sources having any profile whatsoever. 

The algorithms for bus-bars as also for perfect arrays subject to magnetic fields, noise, and/or current drives, all 
involve Green's functions specific to the lattice in question. The Green's function of interest can straightly-forwardly be 
constructed once the eigen-basis of the corresponding connectivity matrix or discrete Laplacian has been determined. 
We have shown that this eigen-basis consists precisely of those linear combinations of eigenfunctions of the connectivity 
matrix for a periodic lattice, which satisfy an appropriate set of boundary conditions. 

The addition and removal of bonds introduces additional boundaries often in the interior of the array. As a result our 
earlier procedure becomes invalid. An island of defects created at the center of the array, e.g., defines a lattice version 
of the sinai billiard problem, which classically has only chaotic solutions. However if the number of bonds eliminated is 
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small, Green's function techniques which have been extensively used to study disordered solids can be very effectively 
applied. We have shown that this techniques yields an 0(N In N + nN) algorithm. The extensions of such algorithms 
to the case of 3D arrays is easy since the corresponding eigen-vectors are specified to be either exp(ik ■ r) or cos^x) 
depending on whether the boundaries are periodic or free boundaries respectively. The eigen-values are now changed 
to 8 — 2cosk x — 2 cos k y — 2 cos k z . In this case, one has to obviously use 3-D transforms. 

It is also worth noting that many of the expressions we have derived can be written in terms of continous functions, 
which are, of course, to be evaluated or sampled only at points belonging to the array. Now if we turn this observation 
around and think of the discrete lattice as being produced by the sampling of the continuous 2-D waveform, we 
make contact with a long-standing problem in electrical engineering, viz. that of the digitisation and subsequent 
recovery of band-width limited analog wave-forms. Usually such waveforms are processed as rectangular spaced 
arrays i. e. they are periodically sampled along the two orthogonally independent variables. However, one can also 
use hexagonal sampling pi] ]. Once discretised, the waveforms can be processed as arrays of numbers x(ni,n 2 ) by 
the computer where (rii,?^) is a discrete point on which the observable x has some value. The periodicity of the 
observable x(n\,n2) is important as it specifies the eigenfunctions associated with the system. The matrix Gq 1 has 
precisely the information about the chosen sampling raster and the periodicity of the lattice. The raster defines 
the total number of nearest neighbours of a general site (711,^2) while the periodicity (or finiteness) of the lattice 
decides the connectivity of the sites at the boundary of the network. The eigenfunctions used in both cases are the 
same. Indeed, the analogy goes deeper and we can think of the JJA itself as a latticized version of the continuum 
case (superconducting islands embedded in a normal background much as the Abrikosov lattice consists of normal 
regions embedded in a superconducting background) . Finally, since hexagonally sampled waveforms can be recovered 
in 0(N In N) steps, this analogy opens up the very interesting possibility of working out O(NlnN) algorithm for a 
triangular JJA. This and other related issues are currently under active investigation. 
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Figure Caption 

Fig.l The figure shows the geometry with currents / being injected and extracted at x — x and x = xl x respectively. 
The • represents the superconducting islands and the line joining them represents a junction. 

Fig. 2 The figure shows the single bus-bar geometry with currents I being injected x = x . The short is at x = Xl x +i- 
Note that one can use any current profile in this case. 

Fig. 3 The figure shows the double bus-bar geometry with shorts at x = xq and x = xl x +i- 
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Fig.l Datta, Das, Sahdev, Mehrotra, Shenoy 




Fig.3 Datta, Das, Sahdev, Mehrotra, Shenoy 



